Convenient Stability Criteria for Difference Approximations of Hyperbolic Initial-Boundary Value Problems

New convenient stability criteria are provided in this paper for a large class of finite-difference approximations to initial-boundary value problems associated with the hyperbolic system u, = Auv + flu + f in the quarter plane x > 0, t > 0. Using the new criteria, stability is easily established for numerous combinations of well-known basic schemes and boundary conditions, thus generalizing many special cases studied in the recent literature. 0. Introduction. In this paper we extend the results of [3] to obtain convenient, more versatile, sufficient stability criteria for a wide class of difference approximations for initial-boundary value problems associated with the hyperbolic system u, = Aux + Bu + f in the quarter plane x > 0, t > 0. Our difference approximation consists of a general difference scheme—explicit or implicit, dissipative or not, two-level or multi-level—and boundary conditions of the type discussed in [3]. We restrict attention to the case where the outflow part of our boundary conditions is translatory; i.e., determined at all boundary points by the same coefficients. Such boundary conditions are commonly used in practice; and in particular, when the boundary consists of a single point, the boundary conditions are translatory by definition. Throughout the paper we assume that the basic scheme is stable for the pure Cauchy problem, and that the assumptions that guarantee the validity of the stability theory of Gustafsson, Kreiss and Sundström [5] hold in our case. With this in mind, we raise the question of stability for our approximation in the sense of Definition 3.3 in [5]. We begin our stability analysis in Secton 1 by recalling Theorem 2.1 of [3], which implies that the entire approximation is stable if and only if the scalar outflow Received September 26, 1983. 1980 Mathematics Subject Classification. Primary 65M10; Secondary 65N10. * Research sponsored in part by the Air Force Office of Scientific Research, Air Force System Command, United States Air Force Grant Nos. AFOSR-79-0127 and AFOSR-83-0150. Part of this research was carried out while the author visited the Institute for Computer Applications in Science and Engineering, NASA Langley Research Center, Hampton, VA 23665. ** Research sponsored in part by the National Aeronautical and Space Administration under NASA Contract No. NASA-17070 while the author was in residence at the Institute for Computer Applications in Science and Engineering, NASA Langley Research Center, Hampton, VA 23665. ©1985 American Mathematical Society 0025-5718/85 $1.00 + $.25 per page 361 License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 362 MOSHE GOLDBERG AND EITAN TADMOR components of its principal part are stable. Thus, our global stability question is reduced to that of a scalar, homogeneous, outflow problem which is the main subject of this paper. Our stability criteria for the reduced problem—stated in the first part of Section 2 and proven in Section 3—depend both on the basic scheme and the boundary conditions, but very little on the intricate interaction between the two. Consequently, our new criteria provide in many cases a convenient alternative to the well-known Gustafsson-Kreiss-Sundstrôm criterion in [5]. In the second part of Section 2, we use the new stability criteria to reestablish all the main examples in our previous paper [3]. We show that if the basic scheme is arbitrary (dissipative or not) and the boundary conditions are generated by either the explicit or implicit right-sided Euler schemes, then overall stability is assured. For dissipative basic schemes we prove stability if the boundary conditions are determined by either oblique extrapolation, by the Box-Scheme, or by the right-sided weighted Euler scheme. Section 2 contains some new examples as well. Among these we find that if the basic scheme is arbitrary and two-level, then horizontal extrapolation at the boundary maintains overall stability. Other stable cases occur when the basic scheme is given by either the Crank-Nicolson scheme or by the backward (implicit) Euler scheme, and the boundary conditions are determined by oblique extrapolation. Such examples, where neither the basic scheme nor the boundary conditions are dissipative, could not have been handled by our previous results in [3]. All told, our examples in this paper incorporate most of the cases discussed in the recent literature; e.g. [l]-[3], [5], [6], [8]-[10], [12], [13], [15], [16]. As in [3], we point out that there is no difficulty in extending our stability criteria to problems with two boundaries. This is so since, if the corresponding left and right quarter-plane problems are stable, then by Theorem 5.4 of [5], the original problem is stable as well. 1. The Difference Approximation and the Reduced Problem. Consider the first-order hyperbolic system of partial differential equations (1.1a) du(x,t)/dt = Adu(x,t)/dx + Bu(x,t)+t(x,t), x > 0, t > 0, where u(x, t) = (m(1)(x, t),. ..,w0, where A1 and Au are of orders / X / and (n I) X (n — I), respectively. The solution of (1.1) is determined uniquely if we prescribe initial values (1.1b) u(x,0) = u°(x), x>0, and boundary conditions (1.1c) uI(0,/) = SuII(0,0 + g(0. *>o, License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use HYPERBOLIC INITIAL-BOUNDARY VALUE PROBLEMS 363 where S is a fixed / X (n — I) matrix, g(/) is a given /-vector, and (1.3) ul = (uw,...,u^Y, un = (W)' is a partition of u into its inflow and outflow components, corresponding to the partition of A in (1.2). In order to solve the initial-boundary value problem (1.1) by a finite-difference approximation we introduce, as usual, a mesh size Ax > 0, At > 0, such that A = At/Ax = constant; and using the notation v„(r) = \(vAx, t), we approximate (1.1a) by a general, consistent, two-sided, solvable, multi-level basic scheme of the form s Q_xyv(t + Ai) = £ Qj,,{t oAt) + Atb,{t), . . £v, = v,+ i' o = -l,...,s,v = r,r +1,.... j—f Here, the n X n matrices A-a are polynomials in A and ArZ?, and the «-vectors b„(r) depend smoothly on f(x, t) and its derivatives. The equations in (1.4a) have a unique solution if we provide initial values (1.4b) y„(oAt) = vv°(oAt), o = 0.s, v = 0,1,2,..., where in addition we must specify, at each time step t = pAt, p = s, s + 1,..., boundary values v„(i + At), v = 0,... ,r 1. As in [3], these boundary values will be determined by two sets of boundary conditions, the first of which is obtained by taking the last n — I components of general boundary conditions of the form T_At + At)= E Tjr.(t oAt) + Atdp(t), 0 = 0 m Ta= HCjaEJ, o--l,...,f,F-0,...,r-l, where the matrices Cja are polynomials in A and ArZ?, the Cj(_X) are nonsingular, and the «-vectors d„(f) are functions of f(x, t) and its derivatives. If we put Cjo = r-n <^iii \ / i \ / ji ciii ciin ' v» Lu ' " in accordance with the partition of A and u in (1.2), (1.3), then this set of conditions takes the explicit form T^Wt + At) + Tnx\u(t + At) (1 4C) = E [Tall\l(t „At) + ro«»v»(r oAt)} + Atd^t), *• ' ' a = 0 m Tjla= E CllaEJ, a = I, II, «v = 0,...,r 1. y=o License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 364 MOSHE GOLDBERG AND EITAN TADMOR We call such boundary conditions translatory since they are determined at all boundary points by the same coefficients. For the second set of boundary conditions we use the analytic condition (1.4d) v0J(r + Ai) = Svll(t + At) + g(t + At), together with r — 1 additional conditions (not necessarily translatory) of the form k \\t + At) = £ [Djp\l(t + At) + Djpnyll(t + At)} + Atel(t), (1.4e) J-1 v = l,...,r1, where the matrices Dj) and Dj11—of orders / X / and / X (n I) respectively—are polynomials in the blocks A" and AtBaß, a, ß = I, II, of the matching partitions Hi ;■)• *HS £)■ so that the D11 vanish whenever B does; and the /-vectors e^i) are functions of f(x, t), g(t), and their derivatives. It has been shown in [3] that Eqs. (1.4c)-(1.4e) can be solved uniquely for the required boundary values v„(t + At), v = 0,.. .,r 1, in terms of neighboring values of v, at least for sufficiently small Ai; and that boundary conditions of the form (1.4e) can be constructed to any degree of accuracy. A concrete example of second-order accurate boundary conditions of the form (1.4e), for the special case B = f = 0, was given in [2]. The difference approximation is completely defined now by (1.4); so assuming that the basic scheme is stable for the pure Cauchy problem (-00 < v < oo), we may ask whether the entire approximation is stable. More precisely, we make from now on the same assumptions about approximation (1.4) as in [5], so that the stability theory of Gustafsson, Kreiss and Sundström holds, and we raise the above stability question in the sense of Definition 3.3 of [5]. As in [3], the first step will be to reduce our stability question to that of a scalar, outflow approximation with homogeneous translatory boundary conditions. This reduction is obtained by applying (1.4a)-(1.4c) to the scalar outflow problem (1.5) u, = aux, a = constant > 0, x > 0, t ^ 0 (which requires no analytic boundary conditions), where a varies over the eigenvalues of Au. In other words, we set A1 = B = i = 0, Au = a = constant > 0; so that (1.1a) yields (1.5); and (1.4a)-(1.4c) reduces to a self-contained, scalar, homogeneous approximation which consists of the basic scheme s Q_xvv(t + At)= £ QMl-oAt), p = r,r+l,..., (1.6a) Qa= E aJOEJ, o = -l,...,s; j—r with initial values (1.6b) vv(oAt) = v°(oAt), a = 0,...,s,v = 0,l,2,...; License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use HYPERBOLIC INITIAL-BOUNDARY VALUE PROBLEMS 365 and translatory outflow boundary conditions i T_xvv(t + At)= £ Tavv(t-oAt), v = 0,...,r-l, (1.6c) m ° = ° Ta=Y,CjaEJ, c0(-d*°o=-l,...,q, 7 = 0 where the scalars aja and cJa are polynomials in a, and the basic scheme (1.6a) is consistent with (1.5). We are now ready to state our main result in Section 2 of [3], which we reformulate as follows: Theorem 1.1 ([3, Theorem 2.1]). Approximation (1.4) is stable if and only if its reduced form in (1.6) is stable for every eigenvalue a of A11. That is, approximation (1.4) is stable if and only if the scalar outflow components of its principal part are stable. The above reduction theorem implies that from now on we may restrict our stability study to the scalar outflow approximation (1.6). Hence, we conclude this section by stating the following five assumptions which will hold throughout the paper, and guarantee the validity of the Gustafsson-Kreiss-Sundström theory [5] for approximation (1.6). Assumption 1.1 ([5, Assumption 3.1]; [3, Assumption 2.1]). Approximation (1.6) is boundedly solvable; i.e., there exists a constant K > 0 such that for each y< e /2(Ax) there is a unique solution w e /2(Ax) to the equations Q-iwv=yr> v = r,r +1,..., T-iK =%> v = 0,...,r 1, with ||h>|| < ZCHj'll; where Q_x and T_x are defined in (1.6a, c), and /2(Ax) is the space of all grid functions w = {w„}f=0 with ||w||2 = AxT^=0\wv\2 < oo. Assumption 1.2 ([5, Assumption 5.1]; [3, Assumption 2.2]). The basic scheme is stable for the pure Cauchy problem -oo < v < oo. That is, if we define the basic characteristic function by (1.7) P(z,k)= £ aj(z)K> j—r where s (1.8) aj(z) = ajl.X) £ z-°-xajo, j = -r,...,p, o = 0 then we have: (i) The von Neumann condition; i.e., the solutions z(k) of the basic characteristic equation (1.9) P(z, k) = 0 satisfy |z(k)| < 1 for all k with \k\ = 1. (ii) If |k| = 1, and if z(k) is a root of (1.9) with |z(k)| = 1, then z(k) is a simple root of (1.9). License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 366 MOSHE GOLDBERG AND EITAN TADMOR Assumption 1.3 (see [5, Assumption 5.4 together with Definition 10.1]; compare [3, Assumption 2.3], and Osher [10]). The basic scheme (1.6a) belongs to the family of schemes for which the Gustafsson-Kreiss-Sundström theory in [5] holds. This family contains in particular the following two classes: (i) Dissipative basic schemes; i.e., schemes for which the roots z(n) of (1.9) satisfy (1.10) |z(k)|<1 for all k with \k\= 1, k + 1. (ii) Unitary basic schemes (also known as strictly nondissipative schemes) where the roots of (1.9) satisfy (1.11) |z(k)| = 1 for all \k\=1. Obviously, if the basic scheme belongs to any of these two classes, then it satisfies the von Neumann condition in Assumption 1.2(i). Assumption 1.4 ([5, Assumption 5.5]; [3, Assumption 2.4]). a_r(z),ap(z) ± 0 forall\z\>l. Assumption 1.5. We assume that m E |c,(z)|#0 forall\z\> 1 where, in analogy with (1.8), (1.12) e/z)^.,,£z —^ j = 0,...,m. o = 0 Assumption 1.5 is necessary for stability, as shown in Remark 3.4 below. This last assumption—which should have been included in [3] as well—is easily verified for all practical boundary conditions. 2. Statement of Results and Examples. In order to state our main stability criteria we define, in complete analogy with (1.7), the boundary characteristic function m R(z,k)= Zcj(z)kJ 7 = 0 where the ciz) are given in (1.12). Defining the function Ü(z,k)=\P(z,k)\ + \R(z, k)\, we shall prove Theorem 2.1 (1st Main Theorem). Approximation (1.6) is stable if (2.1) fi(z,K)*0 for all |z|> 1, 0 < |k| < 1, (z, k) # (1,1). Next, let us divide the (z, k) domain in (2.1) into three disjoint parts and restate Theorem 2.1 as follows: Theorem 2.1' (1st Main Theorem revisited). Approximation (1.6) is stable if (2.2a) ß(z,ic)#0 for all \z\ > 1, |«c| = 1,k# 1; (2.2b) Q(z, k = 1)#0 forall\z\>\,z + \; (2.2c) Q(z,k)*0 for all \z\> 1,0 < |k| < 1. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use HYPERBOLIC INITIAL-BOUNDARY VALUE PROBLEMS 367 The advantage of this new setting is explained by Theorem 2.2 in which we provide useful sufficient conditions for each of the three inequalities in (2.2) to hold. Before stating this theorem, we need the following definitions: Definition 2.1. The boundary conditions (1.6c) are said to be dissipative if the roots z(k) of the boundary characteristic equation (2.3) R(z, k) = 0 satisfy |z(k)|< 1 for a//|/v|< 1,k* 1. Definition 2.2. We say that the boundary conditions (1.6c) satisfy the von Neumann condition if the roots z(k) of (2.3) satisfy |z(«)|< 1 forall\ic\ = 1. Definition 2.3. The boundary conditions (1.6c) are called boundedly solvable if there exists a constant K > 0 so that for each y e /2(Ax) there is a unique solution w e /2(Ax) to (2.4) T_xwr=yv, v = 0,1,2,..., with||w|| 0, v = 0,...,r 1. As in the previous example, the boundary conditions are two-level and dissipative (e.g. [3, Example 3.6]), so Theorem 2.2(i)(iic) gives (2.2a, b). The dissipativity of the boundary conditions also implies the von Neumann condition. And it is trivially verified ([3, Example 3.6] again) that T_x(k) =t 0 for |k| < 1 (where T_x(k) is defined in (2.6)), so that Lemma 2.1(i) implies solvability. Thus, (2.2c) follows from Theorem 2.2(iii); and Theorem 2.T assures stability. Example 2.3 (compare the special cases in [8, Theorem 6], [10, Section XXIII], [5, Theorems 6.1 and 6.3], [1], [6, Theorem 2.1], and [3, Example 3.1]). Take an arbitrary two-level basic scheme, and define the boundary conditions by horizontal extrapolation of order w — 1 : (2.7) vv(t + At)=t [")(-l)J+lv,+J(t + At), v = 0,...,r-l. 7 = 1 W ' We have R(z,k) = R(k) = 1£(7W+V-(1-«)-, 7 = 1 V ' so R(k) ¥= 0 for k =£ 1, which directly gives (2.2a, c). Moreover, since the basic scheme is two-level, Theorem 2.2(iia) implies (2.2b); and Theorem 2.1' again proves stability. It is interesting to note that the above result may fail, both for nondissipative and dissipative basic schemes, if the basic scheme is of more than two time-levels. A nondissipative example was given in Theorem 6.2 of [5] by Gustafsson et al. who License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use HYPERBOLIC INITIAL-BOUNDARY VALUE PROBLEMS 369 showed that the unitary, 3-level Leap-Frog scheme (2.8) vp(t + At) = vp(t-At) + Xa[vp+x(t)-vp„x(t)], v = 1,2,3,..., provides an unstable approximation in combination with the linear boundary extrapolation ((2.7) with u> = 2): v0(t + At) = 2vx(t,+ At) v2(t + At). We proved instability ([3, (3.6)]) for the case where the basic scheme is the 3-level, 5-point, dissipative version of (2.8): vv(t + At) = [z ±(E I)\l ZT1)2]^/ At) + Xa(E E'l)vp(t), 0 < e < 1,0 < Xa < 1 e, v = 2,3,..., and the boundary conditions are given by (2.7) with p = 0,1. Example 2.4 ([3, Example 3.2]; compare the special cases in [2, Example 1] and [6, Theorem 2.2]). Let the basic scheme be dissipative, and determine the boundary conditions by oblique extrapolation of order w — 1: (2.9) vv(t + At) = £ (;)(-l)y+W "O I)*'], " 0.->' ISince the basic scheme is dissipative, then Theorem 2.2(i) implies (2.2a). Furthermore, the boundary characteristic function for (2.9) is r(z,k) = i£(")(-irvv=(i-z-yf, hence (2.10) ß(z, k)>\R(z,k)\*0, z + k. This yields (2.2b, c); and Theorem 2.1' implies stability. Example 2.5 ([3, Example 3.3]; compare the special cases [5, Theorem 6.1], [12, (3.4)], and [2, Example 4]). Take any dissipative basic scheme, and let the boundary conditions be generated by the second-order accurate Box-Scheme v„(t + At) + v„ + x(t + At) Xa[vp+X(t + At) vr(t + At)] = vp(t) + vp+l(t) + Xa[vp + X(t) v„(t)], v = 0,...,r1. As in the previous example, dissipativity implies (2.2a); and Theorem 2.2(iic) implies (2.2b). Next, we recall Example 3.3 of [3] where it was shown without difficulty that the roots z(k) of the boundary characteristic function satisfy |z(«)|= 1 for |k| = 1, and that T_x(k) ¥= 0 for |k| < 1. Thus, the boundary conditions satisfy the von Neumann condition; and Lemma 2.1(i) implies solvability. With this, Theorem 2.2(iii) gives (2.2c), and by Theorem 2.1' stability follows. Example 2.6 ([3, Example 3.4]; compare the special case in [9, Theorem 18.1]). Take any dissipative basic scheme, and define the boundary conditions by the right-sided, 3-level, weighted Euler scheme vr(t + Ai) = vp(t At) + Xa[2vp+X(t) vp(t + At) vp(t At)}, *■ ' 0 < Aa< l,v = 0,...,r1. As in the two previous examples, Theorem 2.2(i) implies (2.2a). Further, we have R(z, k) = 1 z"2 Xa(2KZ'x 1 z~2); License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 370 MOSHE GOLDBERG AND EITAN TADMOR